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The Spalart-Allmaras and the Menter SST k-oj turbulence models are shown to have 
the undesirable characteristic that, for fully turbulent computations, a transition region 
can occur whose extent varies with grid density. Extremely fine two-dimensional grids 
over the front portion of an airfoil are used to demonstrate the effect. As the grid density 
is increased, the laminar region near the nose becomes larger. In the Spalart-Allmaras 
model this behavior is due to convergence to a laminar-behavior fixed point that occurs in 
practice when freestream turbulence is below some threshold. It is the result of a feature 
purposefully added to the original model in conjunction with a special trip function. This 
degenerate fixed point can also cause non-uniqueness regarding where transition initiates 
on a given grid. Consistent fully turbulent results can easily be achieved by either using 
a higher freestream turbulence level or by making a simple change to one of the model 
constants. Two-equation k-ui models, including the SST model, exhibit strong sensitivity 
to numerical resolution near the area where turbulence initiates. Thus, inconsistent ap- 
parent transition behavior with grid refinement in this case does not appear to stem from 
the presence of a degenerate fixed point. Rather, it is a fundamental property of the k-uo 
model itself, and is not easily remedied. 

I. Introduction 

Within the aerodynamics community, the Spalart-Allmaras (SA) 1 and the Menter shear stress 
transport (SST) 2 k-u: turbulence models have been widely-used and trusted for Reynolds-averaged 
Navier-Stokes (RANS) computations for over a decade. The SA model is a one-equation model 
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based primarily on empiricism and on dimensional analysis arguments. It works well for both 
attached and for many separated flows. It was one of the first field-equation models to gain wide 
acceptance into modern aerodynamic application CFD codes. The SST model combined the ro- 
bustness of the k-u turbulence model 3 near walls with the capabilities of the k-e model 4 away from 
walls (also removing the sensitivity to variations in farfield u boundary conditions 5 ). It also in- 
cluded the effects of turbulent shear stress transport through Bradshaw’s assumption 6 that the shear 
stress in a boundary layer is proportional to the turbulent kinetic energy. This “SST” modification 
was the key that gave the model its ability to predict separated flows much better than the standard 
k-u form. Both the SA and SST models have been remarkably successful and stable (very few 
changes proposed or required) over the years. 

Recently, Rumsey et al. 7 showed that certain forms of the low Reynolds number two-equation 
k-e model could exhibit apparent transition behavior that is dependent on initial conditions and 
method of solution (such as different procedures of mesh sequencing). They determined the nu- 
merical reasons behind the long-recognized inconsistent transition prediction behavior of the k-e 
model (see, e.g., Abid 8 ). The apparent laminar flow region predicted by the model upstream of 
transition was termed “pseudo-laminar” because, although the eddy viscosity was predicted to be 
near- zero (yielding laminar results for all intents and purposes), the behavior of the near-zero tur- 
bulent quantities were incorrect relative to each other in the laminar limit. The pseudo-laminar 
state could be reached in some regions of the flow solely because k — e — 0 was a stable solution 
to the equations under certain conditions, and not because of any physics built into the model. In 
a sequel paper, Anders et al. 9 analyzed a broader class of two-equation models, and included the 
effect of a time- varying mean shear. 

CFD practitioners are generally aware that most transport turbulence models “transition” on 
their own. 10-1 2 The location of this apparent transition generally has no relationship with exper- 
imental results. In other words, in computed results there is usually a region of laminar flow 
followed by turbulence, even if a CFD code is run in fully-turbulent mode, with the turbulence 
equations active everywhere in the flow field. For the SA and SST models at reasonably high 
Reynolds numbers, this apparent transition location tends to occur very far upstream in most ap- 
plications, and the models seem to avoid the aforementioned problems often seen with the k-e 
model. 

However, in light of the recent analyses of two-equation models, 7 - 9 it is appropriate to apply 
a similar analysis to the one-equation SA model to determine its ability to produce undesirable 
apparent transition behavior. Also, although we already know from Anders et al. 9 that the standard 
form of the k-u model does not exhibit anomalous behavior (in the sense of converging to a pseudo- 
laminar solution), it is still useful to scrutinize the SST k-u model for transition sensitivities to grid 
density. The main question to try to answer is whether there are circumstances for which the SA 
and SST models can yield inconsistent results in terms of the extent of laminar flow predicted. If 
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present, such an inconsistency can have an impact on the reliability of computed results. 


II. Analysis of the Turbulence Models 

A. Spalart-Allmaras 

The one-equation SA model is written in terms of the turbulence quantity v. 
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The eddy viscosity is given by v t = vf v i, 0 is the magnitude of the vorticity, d is the distance to 
the nearest wall, and the constants are: C bl = 0.1355, C b2 = 0.622, k = 0.41, a = 2/3, C t 3 = 1.2, 
C t4 = 0.5, C w 2 = 0.3, C w 3 = 2, C vl = 7.1, and C',,i = a,/* 2 + (1 + C w )/(7. 

As noted in Spalart and Allmaras, 1 the / /2 term was an ad hoc numerical fix that makes z> = 0 
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a stable solution to the equations with a small basin of attraction (also note that its constants C t 3 
and C'ia changed values between the AIAA paper reference and the subsequent journal article 
reference). The analysis below confirms this characteristic. The reason the model developers 
originally wanted this behavior was due to their use of a “trip function,” f t iAU, not described 
above but given in the original reference. Because the SA model without f t2 would often trip 
to turbulence upstream of where the numerical trip was set for transitional flows, Spalart and 
Allmaras invented the f t2 term to delay it. However, arguably most present-day users of SA do not 
employ the trip function, but rather run the model in fully-turbulent mode (or else force desired 
laminar regions by zeroing out the production term). In the Results section, potential undesirable 
consequences of the presence of the f t 2 term for fully-turbulent computations will be shown. 

For the first part of the analysis here, the homogeneous form of the one-equation model is used. 
As discussed in Rumsey et al ., 7 this is a useful approximation because the log-layer or equilibrium 
layer of a turbulent boundary layer has the same dynamical characteristics as a homogeneous flow. 
The dynamic similarity is exploited by treating the flow as “locally homogeneous.” In other words, 
although the mean shear varies throughout a boundary layer, in the analysis the assumption is made 
of a locally fixed mean shear at a given location. 

The homogeneous form of the SA equation is 
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By inspection, it is clearly seen that v = 0 is a critical point that satisfies the equation in the 
steady-state. The stability of this solution is investigated by solving for dR/dv in the vicinity of 
v = 0 : 
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Using the chain rule and the fact that near v = 0: S — > U, f rl 0, f v2 1, f t2 — »• C t 3 , 
r — > 0 , g — > 0 , f w — > 0 , it can be shown that: 
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Because C t 3 = 1.2, Eq. (13) is always negative, which means that the null critical point is always 
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stable. Therefore, whenever u gets in the vicinity of v = 0, the solution to the partial differential 
equation may be attracted toward D = 0. 

This stability can be demonstrated numerically by solving for the instantaneous time rate-of- 
change of the turbulence quantity for a wide array of possible current state values. This is done by 
first writing Eqs. (10) - (11) in nondimensional form: 
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(The quantity r also requires a (M/Re) factor: r = ( M / Re ) id / ( S' n 2 d' 2 ) . ) In these expressions, 
id = z>/z/ ref , t' = fa re f/L ref , d' = d/L re f , and Q' = QL ve{ /a iei (with a ref representing the reference 
speed of sound). 

Next, a variety of values are input for the independent variables id, d 1 Re/M, and O'. An 
example solution map is shown in Fig. 1. Here, setting fT = 1000, a contour plot of du' /dt' is 
made for an array of u' and d' yj Re/M values. The solution map shows that did /dt' < 0 for the 
region below and to the left of the neutral line (the neutral line is the line along which did / dt' = 0). 
This means that for an initial condition in this region the solution variable v' is driven lower as the 
numerical iteration procedure progresses. Above and to the right of the neutral line, dv' /dt! > 0, 
which means that v' is driven higher. Thus, for d! \J Re/ M to the right of the dot in the figure, 
there are two possible stable solutions. The portion of the neutral line above the dot represents the 
turbulent solution to the homogeneous form of the SA equation. For a given d' Re/M > 0.09 
or so, the solution to the equation will converge to the corresponding u' on the stable neutral line, 
provided that the initial condition is somewhere above the unstable neutral line. Otherwise, the 
solution will converge to the degenerate laminar u' = 0 value. 

Although not shown, for larger d'^j Re/M, the lower part of the neutral line approaches a 
constant value near fd « 0.6, and the upper part of the neutral line continues upward and to the 
right. Also, for different values of O', the general character of the solution map remains the same 
as that shown in the figure, but the specific contour shapes and levels shift. 

This exercise using the homogeneous form of the SA equation is instructive, but not realistic. 
In the full form of the equation, advection and diffusion are locally active and play a significant 
role. These terms can have basically three effects on the solution map, when added to the right- 
hand-side of Eq. (14). First, if the sum of the terms is negative, then the general character of the 
solution map remains unchanged from the homogeneous map: only its levels and location of the 
neutral lines are different. Second, if the sum of the terms is positive and relatively large compared 
to the other terms in the equation, then the possible problem of laminar solution behavior goes 
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away. An example solution map is shown in Fig. 2. For this map, O' is again taken as 1000, but 
this time an assumed nondimensional advection plus diffusion contribution is taken to be a constant 
value of 10, for demonstration purposes. The map shows that in this case there is only one neutral 
line; it is stable and corresponds to the turbulent solution. Any initial condition will converge to 
this turbulent solution. For fully turbulent flow, this situation is the desired one: the numerics are 
such that the solution always converges to the intended, unique value. 

Third, very interesting non-unique behavior occurs if the sum of the right-hand-side advection 
plus diffusion is positive but not too large compared to the other terms. An example solution map 
is shown in Fig. 3. For demonstration purposes, O' is again taken as 1000, and nondimensional 
advection plus diffusion is taken to be a constant value of 3. In this case, there are three distinct 
neutral lines. The lower line oriented horizontally is stable, and will be reached for any initial 
condition below the unstable neutral line. This lower stable neutral line has very low levels of z>', 
yielding even smaller levels of eddy viscosity. Thus, convergence to this line represents a solution 
that behaves like laminar flow. It will be shown in the Results section that this solution behavior 
can cause the apparent transition location for a fully turbulent computation to vary noticeably with 
simple grid refinement. 

B. Menter SST k-u 

The two-equation SST model is written as follows: 
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and the parameters 7, a k , 07, and 6 are calculated via: 
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and 71 = 0.55317, 72 = 0.44035, a kl = 1.17647, a k2 = 1, a ul = 2, a u 2 = 1.16822, ft = 0.075, 
ft = 0.0828, k = 0.41, a L = 0.31, (3* = 0.09, (2 is the magnitude of the vorticity, and d is the 
distance to the nearest wall. Because many CFD codes use the approximation that P k = nSF (as 
recommended by Menter 13 ), it is employed here as well. 

As an initial step, we look briefly at the homogeneous form of the SST model. Because we 
are concerned primarily with transition and the inner part of the boundary layer, we perform this 
initial analysis with F\ = F 2 = 1. Therefore, the cross-diffusion term in Eq. (17) is zero. Nondi- 
mensionalizing the turbulence quantities via: 
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and then, for convenience, further normalizing these variables, along with time, via: 
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the homogeneous form of SST inside the boundary layer can be written: 
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Similar to the SA model, a degenerate solution satisfies these equations in the steady state: 
k* = 0, tv* = This value of u* is less than l/a\ for the values of the given constants, so 

a\k* < k* / uj * , and the stability of the degenerate solution can be found from the eigenvalues of 
the coefficient matrix of the linear system (Eqs. (34) and (35)): 
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Plugging in the degenerate solution, the eigenvalues are found to be real with one positive and one 
negative: this indicates that the solution is a saddle point. 7 - 9 Generally, a saddle point is not stable 
in practice because orbits approach the critical point along one eigenvector, but then recede along 
the eigenvector associated with the unstable solution. Thus, this analysis suggests that the SST 
model should not be attracted to the degenerate solution. 

However, as with the SA model, it is important to also consider the effects of the advective and 
diffusive terms in the SST model in order to truly begin to understand its dynamical behavior. The 
easiest way to do this is to perform an actual computation, and monitor the values of all the terms 
as the solution progresses. This has been done, and detailed results will be shown in the Results 
section. However, for the purposes of analysis and attempting to get a picture of the overall k*- 
u* solution space, we make the assumption that the dominant turbulent transport effects act over 
a distance that represents the length scale characterizing the large turbulent eddies, 14 given by 
£ = k 3 / 2 /e = k 1 / 2 /uj. Then, the nondimensional advection and diffusion terms in the turbulent 
region can be approximated by: 
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w-eqn cross diffusion : 2(1 - ~ 0 {j^J ~ 0 f^' 2 ) = ( C u,,cdWI 41) 

where /i) = /./,/ // re r, and CV/t, C^d, C^d, and C u ,cd are (unknown) coefficients that 

can be determined at instantaneous locations and times in a particular flow field via numerical 
computations. Also, f/j represents the (nondimensional) mean flow velocity components, which 
are taken as locally constant and are absorbed into the coefficients C'k.A and C UtA . Employing the 
normalizations from Eqs. (31) - (33), the inhomogeneous equations can be approximated as: 
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(with F\ and Fo no longer assumed to be 1). These equations are admittedly very crude, because of 
the approximations made in Eqs. (37) - (41). In fact, analysis of the actual advection and diffusion 
terms from CFD results at several locations in the boundary layer confirmed that the coefficients 
C kjA , etc., do not come out to be constant. However, the length scale £ = k 3 ^ 2 /s yielded a fairly 
reasonable approximation that was much better than other choices. The solution to Eqs. (42) and 
(43) can be found. It is the intersection of the A *-nullclinc, given by: 


( i f 

1 - (C k , A )k*-V 2 - C k , D J 


' " F 2 ((3* - (C M )fe‘-V2 - C kiD ) 

and the cc*-nullcline, given by: 


F-2 

u > — 
a\ 

* - F 2 
U < — 
(l\ 


(44) 

(45) 


* = / 7 

^ \P~ (C^, A )k*~^ - C U , D - C UtCD 

The analytic form of the eigenvalues associated with this solution is very complicated, and does not 
yield a clear conclusion regarding the stability of the solution. However, it is easy to numerically 
find the eigenvalues for a point in the turbulent boundary layer of a typical computation. In all 
cases investigated to date, the solution has turned out to be a saddle point. An example phase- 
plane portrait is shown in Fig. 4. This figure shows the k*- and cA-nullclincs using coefficients 
obtained at a point in a converged turbulent boundary layer where nondimensional eddy viscosity 
/i( = 26.47524. In this flow field, Re — lx 10 7 and M = 0.3. The converged solution at this 
point (circled in the figure) has a value of k! = 0.7312424 x 10 -3 , u' = 0.2386361 x 10 -4 , and 
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nondimensional vorticity magnitude O' = 265.3249. Thus k* = 0.1940168 and c o* = 2.998036. 
Other coefficients obtained from the numerical solution were: Ck,n = —0.3368563 x 10“ 2 , C^a = 
-0.4419278 x 10" 2 , C W:D + C^ C d = 0.9705072 x 10“ 2 , C U , A = 0.1760321 x 1(T 2 , F, = 
0.9879624, and F 2 = 0.9999944. The figure also shows the local trajectories (defined by dk* /dt* 
and du* /dt*) of the turbulence variables in phase-space. 

It is clear from the figure that the solution is a saddle point: solution trajectories approach from 
above and below, but recede in both directions along the cu*-nullcline. This is somewhat surprising, 
because at a saddle point an orbit is usually attracted at first but repelled later on. This would seem 
to imply that the SST model is not strictly stable. However, nonlinear and nonlocal effects were 
not accounted for in the analysis, and it turns out that these act to maintain stability at the solution 
point. For example, using a numerical experiment from the example in Fig. 4, when k* at the 
solution point was perturbed slightly to the left, then the local spatial k* derivatives also changed, 
yielding different local advective and diffusive terms which temporarily shifted the intersection 
point of the nullclines even further to the left; as a result the solution was driven back to the right, 
toward its original position. In any case, SST certainly is numerically stable in practice. 

Although not shown here, the standard k-co model 3 yields a phase-plane portrait very similar to 
that of SST, with its solution also appearing as a saddle point. Thus, its behavior is expected to be 
similar. 

Now the following question arises: what are the effects of this model’s seemingly tenuous 
linearized- solution stability on its apparent transition behavior? The answer is not entirely clear 
from this analysis, but in the Results section some details will be shown that will hopefully shed 
some light on this question. 


III. Numerical Method 

The computer code CFL3D 15 solves the three-dimensional, time-dependent compressible RANS 
equations with an upwind finite-volume formulation (it can also be exercised in two-dimensional 
mode of operation for 2-D cases). Upwind-biased third-order spatial differencing is used for the 
inviscid terms, and viscous terms are centrally differenced. The code originally solved the thin- 
layer form of the equations (in each coordinate direction), but the full Navier-Stokes terms (i.e., 
cross-derivative terms) have recently been added. All solutions shown below use the full Navier- 
Stokes terms, although thin-layer has also been employed and makes almost no difference for the 
particular cases run. 

The CFL3D code is advanced in time with an implicit approximate factorization method. The 
implicit derivatives are written as spatially first-order accurate, which results in block tridiagonal 
inversions for each sweep. However, for solutions that utilize Roe flux-difference splitting, 16 the 
block tridiagonal inversions are further simplified using a diagonal algorithm with a spectral radius 
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scaling of the viscous terms. 

The turbulence models, including SA and SST, are solved uncoupled from the mean flow equa- 
tions using implicit approximate factorization. Their advective terms can be solved using either 
first-order or second-order upwind differencing, with first-order the default for the code. 

IV. Results 

In order to demonstrate inconsistent apparent transition behavior of the SA and SST models, 
it was necessary to perform computations on very fine grids. Easily noticeable inconsistencies did 
not show up using coarser grids. The model problem chosen was the computation over the front 
upper 45 % of a NACA 0012 airfoil, at M = 0.3, Re = 10' based on total chord length, and 
a = 0°. Fully-turbulent computations were performed, which means that the turbulence model 
equations were active everywhere in the flow field. The freestream eddy viscosity was set using a 
typical value 17 of turbulence Reynolds number Rt,oo = ^/(^oo^oo) = 0.1, yielding /j, t>00 = 0.009. 
For the SA model, this was achieved using = 1.341946. It will be shown below that for the 
SA model, the freestream turbulence level chosen can be an important consideration. For the SST 
model, the freestream turbulent quantities used were: k' OQ = 9 x 10 -9 and c o'^ = 1 x HP 6 . 

The finest grid was 769 x 1025, with coarser grids created by successively removing every 
other grid point from the next finest level. Thus the coarser grid sizes were: 385 x 513, 193 x 257, 
97 x 129, and 49 x 65. The grids had a farfield extent of 50 chords, and the minimum normal spacing 
at the wall yielded an average y+ level of 0.05 for the 769 x 1025 grid, 0.1 for the 385 x 513 grid, 
0.2 for the 193 x 257 grid, 0.4 for the 97 x 129 grid, and 0.8 for the 49 x 65 grid. 

A picture of the coarsest grid is shown in Fig. 5 . Note that this coarsest level is on a par with 
typical grid resolution used for many turbulent CFD computations over airfoils and wings. For 
boundary conditions, the part of the grid coming into the nose of the airfoil was set as a symme- 
try plane, the airfoil body used no-slip adiabatic wall conditions, the downstream exit plane used 
extrapolation, and the farfield (upstream and above the airfoil) used a farfield Riemann invariant 
condition. For most of the computations, each grid level was run by restarting from the previous 
coarse grid solution, as is typically done for full- approximation scheme (FAS) multigrid simula- 
tions. 18 All results were converged 5-6 orders of magnitude on each grid (the final density residual 
dropped to 10 _1 ° - 1CT 11 ). 

A. Spalart-Allmaras 

Fig. 6 shows computed surface skin friction coefficients on the airfoil using SA for each of the 
grids. As the grid was refined, each successive solution produced a larger “laminar region” prior 
to going turbulent. Only on the finest 769 x 1025 grid did the transition prediction appear to stop 
changing significantly. For the coarsest grid, the nondimensional eddy viscosity first exceeded 1.0 
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in the boundary layer near x/c = 0.0044. For successively finer grids, the location was 0.0056, 
0.0099, 0.0111, and 0.0112, respectively. As will be shown next, this apparent later transition on 
the finer grids is a numerical artifact related to the existence of the laminar-behavior stable neutral 
line for this turbulence model. 

Solution maps can be created for this realistic computation by monitoring the various in- 
dependent quantities during the course of the solution’s progress. In this case, the quantities 
were recorded during the solution on the 193 x 257 grid, at a point in the boundary layer near 
x/c = 0.00694, where d 1 y Re/M = 0.1157. At this position, the converged solution on the 
97 x 129 grid was turbulent (i.e., the peak nondimensional eddy viscosity level - located some- 
what further from the wall at this location - exceeded 1), but during iterations on the 193 x 257 
grid this location became laminar. 

Fig. 7 shows the descent of the V quantity at this location from near 3.9 initial value to a con- 
verged value of 0.205. (These correspond to nondimensional eddy viscosity levels of about 0.48 
and 4 x 10“ 6 , respectively). This solution on the 193 x 257 was restarted from the converged 
solution on the 97 x 129 grid at multigrid cycle number 8001. Next, by using the actual instan- 
taneous values for O' and nondimensional advection plus diffusion at this location, the solution 
maps were plotted at multigrid cycle numbers 8940, 9120, and 1 1000 (at cycle 8940: O' = 3.1064, 
O' = 5365.3, and advection plus diffusion = —748.8; at cycle 9120: 0' = 2.313, O' = 5399.4, and 
advection plus diffusion = —1137.4; at cycle 11000: O' = 0.205, O' = 5500.3, and advection plus 
diffusion = 24.5). These maps are shown in Figs. 8, 9, and 10. At cycle 8940, the solution map 
shows that the location of the current value of O' (indicated by the open circle) lies slightly above 
and to the left of the stable neutral line in the map space. Thus, its value is being driven lower with 
successive iterations. As its value decreases, the local advection and diffusion also change signifi- 
cantly such that the stable neutral line is driven further to the right (Fig. 9). Finally, at convergence 
(Fig. 10), the final value of O' lies on the stable neutral line representative of a laminar-behavior 
solution. 

Because two stable attractors exist, it is possible to obtain different solutions with SA depend- 
ing on numerical factors such as initial condition (I.C.) and method of solution. An example is 
shown for the 97 x 129 grid in Fig. 11. The result labeled “I.C. #1” is the same result on this grid 
level shown earlier in Fig. 6. The result labeled “I.C. #2” was obtained from an initial condition 
using a converged result on the 193 x 257 grid interpolated onto the coarser level, then re-iterated 
to convergence. Two different apparent transition locations result: the first near x/c = 0.0056 and 
the second near x/c = 0.0097. 

The numerically-induced laminar behavior of the SA model is easily avoided in two different 
ways. The first method, as suggested by Spalart, 11 is to employ higher levels of freestream turbu- 
lence. For example, using z>^ = 3 (/4 = 0.2104), the problem of inconsistent transition was 

eliminated, as shown in Fig. 12. The apparent transition location remained near the airfoil leading 
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edge (near x/c = 0.004 at this Reynolds number) for all grid levels. 

The earlier SA analysis also suggests a second method for avoiding the problem. By setting 
the constant C t 3 = 0, the right-hand side of Eq. (13) remains always positive, and the null critical 
point is unstable (z> = 0 is not an attractor). The same NACA 0012 airfoil case was run with 
n' t oo = 0.009 and C t3 = 0. Although not shown, results again show consistent behavior, and skin 
friction results appear identical to those shown in Fig. 12. 

B. Menter SST k-uj 

The same NACA 0012 case was run using the SST model in fully turbulent mode. Results are 
shown in Fig. 13, and are similar to those using SA: as the grid was refined, the apparent transition 
location moved further aft. Even on the finest 769 x 1025 grid, the location still appeared to 
be influenced by the grid density. For the coarsest grid, the nondimensional eddy viscosity first 
exceeded 1.0 in the boundary layer near x/c = 0.0054. For successively finer grids, the location 
was 0.0062, 0.0076, 0.0086, and 0.0089, respectively. It should be noted that, unlike the SA model 
above and the k-s model reported in Rumsey et al., 7 SST did not appear to exhibit sensitivity to 
initial conditions or method of solution. In other words, a given grid always gave the same result. 

The dynamics of the SST solution behavior near the leading edge region was investigated by 
monitoring various quantities during the solution progress on the 193 x 257 grid, at a point in the 
boundary layer near x/c = 0.00694 (distance from the wall was d!^J Re/ M = 0.2285). At this 
streamwise location, the solution on the 97 x 129 grid was turbulent (peak nondimensional eddy 
viscosity exceeded 1 in the boundary layer), but during iterations on the 193 x 257 grid the peak 
eddy viscosity dropped below 1. 

Fig. 14 shows the changes in n[, k*, and uj* through iteration to convergence. The value of u* 
changed relatively little, but k* decreased by nearly 75%, causing a similar dramatic decrease in 
eddy viscosity. At multigrid cycle number 8500, the phase-plane portrait in Fig. 15 shows that the 
unconverged values of k* and u* (indicated by the open circle) lie to the left of the intersection 
point of the nullclines. Thus, the solution in phase-space continues to be driven to the left. As it 
moves, the nullclines also change, but their intersection remains to the right of the solution until 
convergence is reached. Finally, at convergence (Fig. 16), the solution and the intersection of the 
nullclines coincide, as they must. 

Unlike SA, the SST model continued to exhibit inconsistent transition results even when the 
freestream turbulence level was raised to levels above 1. For example, Fig. 17 shows surface skin 
friction coefficients with fi' too = 1.294 (k'^ = 1.294 x 10 -6 ). Compared with results in Fig. 13, 
these results using higher freestream turbulence show less grid influence, but the influence is still 
visible nonetheless. 

From these examples, it seems evident that the SST model does indeed exhibit a sensitivity of 
its predicted transition location to grid density, but this sensitivity is not due to a null critical point 
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acting as an additional attractor. Rather, there is only one possible solution, and the model dynam- 
ical behavior itself (near the area where turbulence initiates) appears to be particularly sensitive 
to small changes induced by different numerical resolution. One of the reasons for this sensitivity 
may be the fact that, in a linear sense, the solution is a saddle point in k*-LU* phase-space. No 
simple cure was found to remedy this problem. 

V. Conclusions 

In conclusion, both the Spalart-Allmaras and the Menter SST k-uj turbulence models were 
shown to have the undesirable characteristic that, for fully turbulent computations, their apparent 
transition location can vary with grid density. As the grid was refined for an example airfoil com- 
putation, an apparent laminar region near the nose of the airfoil grew larger. In the SA model this 
behavior was due to a feature purposefully added to the original model in conjunction with a spe- 
cial trip function, which many people employing the model do not even use. Solution maps were 
used to show that this feature can cause convergence to a laminar-behaving fixed point. Because 
two stable fixed points exist, even a given grid can produce non-unique results regarding where 
transition initiates. It was also shown that consistent fully turbulent SA results can be achieved in 
practice by either using higher freestream turbulence levels (z>^ = 3 or greater), or by setting the 
constant C t 3 in the model to 0. 

The SST model’s dynamical behavior (near the area where turbulence initiates) appears to 
be particularly sensitive to small changes induced by changes in numerical resolution. Thus, its 
apparent transition location change with grid refinement seems to be a fundamental property of the 
k-oj model itself, and is not easily remedied. Raising freestream turbulence levels improved but 
did not fix the problem. Phase-plane portraits, which included the nullclines of the k-u equations, 
were used to demonstrate the SST solution behavior. This model does not suffer from the problem 
of possible non-uniqueness on a given grid. 

It is important to note that the problems illustrated in this paper do not appear to be very sig- 
nificant at high Reynolds numbers typically employed for aerodynamic calculations. For example, 
the transition location movement at Re = 10 million was less than 1% chord for the NACA 0012 
airfoil. In other words, generally these two turbulence models tend to transition near enough to the 
leading edge of airfoils and wings so that the overall solution remains globally consistent. How- 
ever, this may be problem and/or Reynolds number dependent. In any case, running the SA model 
with z/^ > 3 for computations that do not use the trip function appears to be a good idea. Setting 
C t 3 = 0 is also a viable choice, but may be less desirable from the point of view of model backward 
compatibility, because it is possible that the model behavior in and near separation could also be 
affected by this change. For k-uj models, including SST, users should be aware of the apparent 
transition-location dependence on grid density. In particular, as computers become more powerful 
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and finer grids are more routinely employed, the problem behavior may manifest itself to a greater 
degree. 
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Figure 1. Homogeneous SA solution map of dv'/dt' contours for a range of initial conditions for v' and 

d'y/Re/M ; Cl' = 1000. 
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Figure 2. Inhomogeneous SA solution map of dv' jd t! contours for a range of initial conditions for v' and 
d' sjRe/M; i Y = 1000; advection plus diffusion contribution taken to be 10. 
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Figure 4. Inhomogeneous phase-plane portrait of the SST model in k* and co* space for a typical turbulent 
solution point (circled). 
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Figure 5. Portion of coarsest NACA 0012 grid, 49 x 65. 
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Figure 7. Convergence of v' at specific location in the boundary layer on 193 x 257 grid, SA model with 
n't co = 0.009; solution maps will be plotted for the three times indicated in this figure. 
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Figure 8. SA solution map for NACA 0012 case on 193 x 257 grid at a specific point in the boundary layer near 
x/c = 0.00694, at multigrid cycle number 8940; open circle indicates the particular (unconverged) solution at 
this time. 
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Figure 9. SA solution map for NACA 0012 case on 193 x 257 grid at a specific point in the boundary layer near 
x/c = 0.00694, at multigrid cycle number 9120; open circle indicates the particular (unconverged) solution at 
this time. 
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Figure 10. SA solution map for NACA 0012 case on 193 x 257 grid at a specific point in the boundary layer near 
x/c = 0.00694, at multigrid cycle number 11000; open circle indicates the particular (converged) final solution. 
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Figure 12. Computed surface skin friction coefficient for the NACA 0012 airfoil on various grids, SA model 
with /4,oo = 0.2104. 
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Figure 13, 

with /4,oo 


Computed surface skin friction coefficient for the NACA 0012 airfoil on various grids, SST model 

= 0.009. 
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Figure 14. Convergence of fi' t , k*, and to* at specific location in the boundary layer on 193 x 257 grid, SST 
model with n' too = 0.009; phase-plane portraits will be plotted for the two times indicated in this figure. 
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Figure 15. SST phase-plane portrait for NACA 0012 case on 193 x 257 grid at a specific point in the boundary 
layer near x/c = 0.00694, at multigrid cycle number 8500; open circle indicates the particular (unconverged) 
solution at this time. 
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Figure 16. SST phase-plane portrait for NACA 0012 case on 193 x 257 grid at a specific point in the boundary 
layer near x/c = 0.00694, at multigrid cycle number 11000; open circle indicates the particular (converged) 
final solution. 
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Figure 17. 

with /4,oo 


Computed surface skin friction coefficient for the NACA 0012 airfoil on various grids, SST model 

= 1.294. 
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